R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(osmdata)

Warning: package ‘osmdata’ was built under R version 4.2.1
Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

library(sf)

Warning: package ‘sf’ was built under R version 4.2.1
Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

library(ggplot2)

Warning: package ‘ggplot2’ was built under R version 4.2.1

library(plotly)

Warning: package ‘plotly’ was built under R version 4.2.1
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(rgeos)

Warning: package ‘rgeos’ was built under R version 4.2.1
Loading required package: sp
Warning: package ‘sp’ was built under R version 4.2.1
rgeos version: 0.5-9, (SVN revision 684)
 GEOS runtime version: 3.9.1-CAPI-1.14.2 
 Please note that rgeos will be retired by the end of 2023,
plan transition to sf functions using GEOS at your earliest convenience.
 GEOS using OverlayNG
 Linking to sp version: 1.5-0 
 Polygon checking: TRUE 

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ tibble  3.1.7     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1
✔ purrr   0.3.4     
── Conflicts ────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks plotly::filter(), stats::filter()
✖ dplyr::lag()    masks stats::lag()

get_overpass_url()

[1] "https://overpass.kumi.systems/api/interpreter"

FYI all osmdata object can be combined using the c() operator. This c() operation can only be performed on simple feature (sf) type geographic features. I.e. you have to use osmdata_sf() to create the simple feature or convert sp type features to sf type features. I think that calling the bbox based on the county perimeter is better. There are some campuses that are not labeled only “University of Georgia” explicity, so it may be best to use st_… methods to get landuse polygons based on the UGA operator key and then use those polygons to filter and select UGA buildings. I could also use the trim_osmdata() function to select from the list of all UGA operated landuse polygons and run through them, one-by-one and obtain a building list that way. I need to use format_out = “polygon” in the getbb() function in order for this to work. There may be some edge-cases (e.g. building is on UGA campus landuse polygon, yet is not a UGA operated building.) that lead me to think that simply adding operator=“University of Georgia” to each building is worth the effort.

Pro-tip: value_exact needs to be false in order for match_case to also be false and not fail the query.

bb_county <- getbb('Athens-Clarke County')
all_Buildings <- opq(bbox = bb_county) %>%
  add_osm_feature(key = "building") %>%
  add_osm_feature(key = "operator", value = "University of Georgia", value_exact = FALSE, match_case = FALSE) %>%
  osmdata_sf()

ugaCampus <- opq(bbox = bb_county) %>%
  add_osm_feature(key = "landuse", value = "property") %>%
  add_osm_feature(key = "name", value = "University of Georgia", value_exact = FALSE, match_case = FALSE) %>%
  osmdata_sf()
#Create sf dataset of the UGA campus polygons.
ugaCampus <- ugaCampus$osm_polygons

This query allows me to access the multipolygon that comprises the UGA campus. I will then use this to create a bbox based on this multipolygon. I need to also filter out only the polygons from this multipolygon that are explicitly labeled as universities (i.e. amenity = university).

ugaCampus <- ugaCampus %>% filter(amenity == "university")

Get all bicycle parking features in Athens. There are, I believe, at least two bicycle parking features mapped as polygons, so in order for a completely accurate bike parking count, I need to include them. It may cause some issues downstream.

all_Bicycle_Nodes <- opq(bbox =bb_county) %>%
  add_osm_feature(key = "amenity", value = "bicycle_parking") %>%
  osmdata_sf()
all_Bicycle_Nodes <- all_Bicycle_Nodes$osm_points

Filter only bicycle nodes within the UGA campus.

mat = st_intersects(all_Bicycle_Nodes, ugaCampus, sparse = FALSE)

mat <- apply(mat, 1, any)

all_Bicycle_Nodes_Trimmed <- all_Bicycle_Nodes[mat,]

Export the bicycles nodes and UGA buildings as shapefiles.

st_write(all_Buildings,
         "UGA_buildings.shp",
         layer_options = "ENCODING=UTF-8",
         delete_layer = TRUE)

Warning in abbreviate_shapefile_names(obj) :
  Field names abbreviated for ESRI Shapefile driver

Deleting layer `UGA_buildings' using driver `ESRI Shapefile'
Writing layer `UGA_buildings' to data source `UGA_buildings.shp' using driver `ESRI Shapefile'
options:        ENCODING=UTF-8 
Writing 437 features with 68 fields and geometry type Polygon.

st_write(all_Bicycle_Nodes_Trimmed,
         "Bicycle_parking.shp",
         layer_options = "ENCODING=UTF-8",
         delete_layer = TRUE)

Warning in abbreviate_shapefile_names(obj) :
  Field names abbreviated for ESRI Shapefile driver

Deleting layer `Bicycle_parking' using driver `ESRI Shapefile'
Writing layer `Bicycle_parking' to data source `Bicycle_parking.shp' using driver `ESRI Shapefile'
options:        ENCODING=UTF-8 
Writing 211 features with 22 fields and geometry type Point.

Combine the UGA bicycle parking and buildings datasets.

building_Bike_Combined_Data <- c(all_Buildings, all_Bicycle_Nodes_Trimmed)

ggplotly interactive map of UGA buildings and bike racks.

q <- ggplot() + geom_sf(data = building_Bike_Combined_Data$osm_points,
  fill = 'light blue') + theme_minimal() + geom_sf(data = building_Bike_Combined_Data$osm_polygons,
  fill = 'light blue') + theme_minimal()
ggplotly()

b is building c is bike nodes. Think code is only able to do simple polygons at the time. I will need to add multi-polygon ability in the future.

b <- building_Bike_Combined_Data$osm_polygons
c <- building_Bike_Combined_Data$osm_points

##  First project data into a planar coordinate system (here EPSG:2240)
utmStr <- "+proj=EPSG:2240 +zone=%d +datum=WGS84 +units=ft +no_defs +ellps=GRS80"
crs <- st_crs(2240)

#Perform spatial transform on bike nodes and building polygons
pUTM <- st_transform(b, crs)
ptsUTM <- st_transform(c, crs)

## Set up containers for results
n <- nrow(ptsUTM)
nearestBuilding <- character(n)
distToNearestBuilding <- numeric(n)
buildingCapacity <- numeric(n)
buildingCount <- numeric(n)
buildingRackType <- character(n)
buildingRackCovered <- character(n)

The st_nearest_feature function may also be really useful here. I think I would like to add a column to the original unprojected simple feature type.

## For each point, find name of nearest polygon (in this case, UGA buildings)
for (i in seq_along(nearestBuilding)) {
    buildingCapacity <- ptsUTM$capacity[i,]
    gDists <- st_distance(ptsUTM[i,], pUTM, byid=TRUE)
    nearestBuilding[i] <- pUTM$name[which.min(gDists)]
    distToNearestBuilding[i] <- min(gDists)
}

## Check that it worked
nearestDistanceName <- data.frame(nearestBuilding, distToNearestBuilding)

#E.g. should return something like this.
#       nearestCanton distToNearestCanton
# 1             Wiltz           15342.222
# 2        Echternach            7470.728
# 3            Remich           20520.800
# 4          Clervaux            6658.167
# 5        Echternach           22177.771
# 6          Clervaux           26388.388
# 7           Redange            8135.764
# 8            Remich            2199.394
# 9  Esch-sur-Alzette           11776.534
# 10           Remich           14998.204

plot(c, pch=16, col="red")
#text(c, 1:10, pos=3)
plot(b, add=TRUE)
text(b, b$name, cex=0.7)

Write the dataframe to a .csv file.

write(nearestDistanceName, "test.dbf")

Error in cat(x, file = file, sep = c(rep.int(sep, ncolumns - 1), "\n"),  : 
  argument 1 (type 'list') cannot be handled by 'cat'

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.